Libraries and paths

library(easypackages)
libraries("here","tidyverse","ggeasy","patchwork","lmerTest","emmeans","reshape2")

plotdir = here("plots")
fontSize = 16

baseline_end = 15
transition_end = 35

# cols2use = c("#619cff","#156cff")
cols2use = c("#619cff","#ffa500")

Functions to use later

# functions to use for plotting
make_timeplot <- function(data4plot, x_var, y_var, grp_var,
                          baseline_end=15, transition_end=35,
                          yLabel,title2use,fontSize=16){
  
  p = ggplot(data = data4plot,
           aes_string(x = x_var, y = y_var, group = grp_var)) +
  facet_grid(. ~ grp) 
  
  p = p + geom_smooth(se = FALSE,
                    colour='gray75',
                    alpha=0.1)
  
  p = p + geom_smooth(se=TRUE, aes(group = interaction(grp,condition),
                                 colour = condition,
                                 fill = condition),
                    size=3)
  
  p = p + geom_vline(xintercept = (baseline_end)) +
    geom_vline(xintercept = transition_end) +
    geom_hline(yintercept = 0) +
    xlab("Time Window") + 
    ylab(sprintf("Baseline Normalized %s",yLabel)) +
    guides(colour="none", fill="none") +
    ggtitle(title2use) + easy_center_title() +
    theme_bw() + 
    theme(plot.title = element_text(hjust = 0.5, size=fontSize),
      strip.text.y = element_text(size=fontSize),
      axis.text.x = element_text(size=fontSize),
      axis.text.y = element_text(size=fontSize-2),
      axis.title.x = element_text(size=fontSize),
      axis.title.y = element_text(size=fontSize),
      strip.text = element_text(size = fontSize))
  
  return(p)
  
} # make_timeplot

# make_scatterplot
make_scatterplot <- function(data4plot, x_var, y_var, grp_var,
                             cols2use,xLabel,yLabel,title2use, dotSize=3){
  p = ggplot(data = data4plot, aes_string(x = x_var, y = y_var, colour=grp_var)) + 
  geom_point(size = dotSize) + 
  geom_smooth(method=lm, se=FALSE) + 
  scale_colour_manual(values=cols2use) + 
  xlab(sprintf("Baseline Normalized %s",xLabel)) +
  ylab(sprintf("Baseline Normalized %s",yLabel)) +
  ggtitle(title2use) + easy_center_title() +
  guides(colour="none") +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, size=fontSize),
    strip.text.y = element_text(size=fontSize),
    axis.text.x = element_text(size=fontSize),
    axis.text.y = element_text(size=fontSize-2),
    axis.title.x = element_text(size=fontSize),
    axis.title.y = element_text(size=fontSize),
    strip.text = element_text(size = fontSize))

  return(p)
} # function make_scatterplot

# make_blankplot
make_blankplot <- function(exp, dtype="H", lineSize = 3,ylimits =c(-6,6),fontSize = 16){
  
  data_FR = data.frame(time = c(0,15,16,35,36,60, 0,15,16,35,36,60),
                    condition = c("baseline","baseline",
                                  "transition","transition",
                                  "treatment","treatment",
                                  "baseline","baseline",
                                  "transition","transition",
                                  "treatment","treatment"),
                    grp = c("DREADD","DREADD","DREADD","DREADD","DREADD","DREADD",
                            "SHAM","SHAM","SHAM","SHAM","SHAM","SHAM"))
  data_H = data_FR
  
  if (dtype=="H"){
    
    if (exp=="camk"){
      data_FR$score = c(0,0,0,3,3,3, 0,0,0,0,0,0)
      data_H$score = c(0,0,0,-3,-3,-3, 0,0,0,0,0,0)
      title2use = "CamkII-hM3D(Gq)"
    }else if (exp=="hm4di"){
      data_FR$score = c(0,0,0,-3,-3,-3, 0,0,0,0,0,0)
      data_H$score =  c(0,0,0,3,3,3, 0,0,0,0,0,0)
      title2use = "hM4Di"
    }
    data_H$facet_var = "H"
    data_FR$facet_var = "Firing Rate"
    
  } else{
    
    if (exp=="camk"){
      data_H$score = c(0,0,0,3,3,3, 0,0,0,0,0,0)
      data_FR$score = c(0,0,0,3,3,3, 0,0,0,0,0,0)
      title2use = "CamkII-hM3D(Gq)"
    }else if (exp=="hm4di"){
      data_H$score = c(0,0,0,-3,-3,-3, 0,0,0,0,0,0)
      data_FR$score =  c(0,0,0,-3,-3,-3, 0,0,0,0,0,0)
      title2use = "hM4Di"
    }
    data_H$facet_var = dtype
    data_FR$facet_var = "Firing Rate"

  }
  
  data = rbind(data_FR,data_H)
  data$facet_var = factor(data$facet_var, levels = c("Firing Rate", dtype))
  data$grp = factor(data$grp, levels = c("DREADD","SHAM"))
  data$condition = factor(data$condition, levels = c("baseline","transition","treatment"))
  
  p = ggplot(data = data, aes(x = time, y = score, colour = condition)) + facet_grid(. ~ facet_var) +
      geom_hline(yintercept = 0) +
    geom_line(aes(linetype=grp),size = lineSize) + 
    scale_linetype_manual(values = c("solid", "dashed")) +
    # xlim(xlimits[1],xlimits[2]) +
    scale_y_continuous(limits = ylimits) +
    xlab("Time Window") + ylab("Baseline Normalized Value") +
    geom_vline(xintercept = 15) +
    geom_vline(xintercept = 35) +
    guides(colour="none", linetype="none") +
    ggtitle(title2use) + easy_center_title() +
    theme_bw() +   
    theme(plot.title = element_text(hjust = 0.5, size=fontSize),
      strip.text.y = element_text(size=fontSize),
      axis.text.x = element_text(size=fontSize),
      axis.text.y = element_text(size=fontSize-2),
      axis.title.x = element_text(size=fontSize),
      axis.title.y = element_text(size=fontSize),
      strip.text = element_text(size = fontSize))
  
  return(p)
}

# run_lme
run_lme <- function(data, dv_var){
  
  # formula
  form2use = as.formula(sprintf("%s ~ grp*condition + (1|id)", dv_var))
  
  # run the model
  mod2use = lmer(formula = form2use, data = data)
  aov_tab = anova(mod2use)

  pairwise_res = emmeans(mod2use, pairwise ~ grp*condition)
  
  result = list(model = mod2use,
                aov_tab = aov_tab, 
                pairwise_res = pairwise_res)
  return(result)
  
} # function run_lme

Theoretical plots

CamkII-hM3D(Gq)

camk_plot = make_blankplot(exp = "camk", dtype = "H", fontSize=24)
ggsave(filename = file.path(plotdir,"camk_H_theoretical_prediction.pdf"), plot = camk_plot)
camk_plot

camk_plot = make_blankplot(exp = "camk", dtype = "1/f slope", fontSize=24)
ggsave(filename = file.path(plotdir,"camk_slope_theoretical_prediction.pdf"), plot = camk_plot)
camk_plot

camk_plot = make_blankplot(exp = "camk", dtype = "Broadband Power", fontSize=24)
ggsave(filename = file.path(plotdir,"camk_power_theoretical_prediction.pdf"), plot = camk_plot)
camk_plot

hM4Di

hm4di_plot = make_blankplot(exp = "hm4di", dtype = "H", fontSize=24)
ggsave(filename = file.path(plotdir,"hm4di_H_theoretical_prediction.pdf"), plot = hm4di_plot)
hm4di_plot

hm4di_plot = make_blankplot(exp = "hm4di", dtype = "1/f slope", fontSize=24)
ggsave(filename = file.path(plotdir,"hm4di_slope_theoretical_prediction.pdf"), plot = hm4di_plot)
hm4di_plot

hm4di_plot = make_blankplot(exp = "hm4di", dtype = "Broadband Power", fontSize=24)
ggsave(filename = file.path(plotdir,"hm4di_power_theoretical_prediction.pdf"), plot = hm4di_plot)
hm4di_plot

CamkII-hM3D(Gq) Firing Rate

data_file = here("experimental_data","camk_databasenorm.csv")
data = read.csv(data_file)
data$grp[data$grp=="exp"] = "DREADD"
data$grp[data$grp=="ctrl"] = "SHAM"
data$grp = factor(data$grp)
data2use = data %>% 
  filter(condition=="treatment")

# Firing Rate
p0 = make_timeplot(data4plot=data, x_var="time", y_var="r_basenorm", grp_var="id",
                   yLabel="Firing Rate",title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_FR_timeplot.pdf"), plot = p0)
p0

# run LME
res = run_lme(data = data, dv_var = "r_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)    
## grp            21.56  21.559     1  11.02  7.3291 0.02036 *  
## condition     247.81 123.906     2 763.00 42.1230 < 2e-16 ***
## grp:condition 378.05 189.026     2 763.00 64.2610 < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 4.440802e-18
res$aov_tab[3,6]
## [1] 1.613408e-26
res$pairwise_res
## $emmeans
##  grp    condition  emmean    SE   df lower.CL upper.CL
##  DREADD baseline    0.000 0.644 12.7    -1.39    1.393
##  SHAM   baseline    0.000 0.509 12.7    -1.10    1.102
##  DREADD transition  2.759 0.636 12.2     1.38    4.143
##  SHAM   transition -0.147 0.503 12.2    -1.24    0.947
##  DREADD treatment   3.122 0.631 11.8     1.74    4.500
##  SHAM   treatment  -0.399 0.499 11.8    -1.49    0.690
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                             estimate    SE    df t.ratio p.value
##  DREADD baseline - SHAM baseline         0.000 0.821  12.7   0.000  1.0000
##  DREADD baseline - DREADD transition    -2.759 0.262 763.0 -10.531  <.0001
##  DREADD baseline - SHAM transition       0.147 0.817  12.5   0.180  1.0000
##  DREADD baseline - DREADD treatment     -3.122 0.251 763.0 -12.462  <.0001
##  DREADD baseline - SHAM treatment        0.399 0.815  12.4   0.490  0.9957
##  SHAM baseline - DREADD transition      -2.759 0.815  12.4  -3.387  0.0463
##  SHAM baseline - SHAM transition         0.147 0.207 763.0   0.710  0.9808
##  SHAM baseline - DREADD treatment       -3.122 0.811  12.2  -3.850  0.0215
##  SHAM baseline - SHAM treatment          0.399 0.198 763.0   2.016  0.3339
##  DREADD transition - SHAM transition     2.906 0.811  12.2   3.584  0.0338
##  DREADD transition - DREADD treatment   -0.363 0.230 763.0  -1.577  0.6142
##  DREADD transition - SHAM treatment      3.158 0.809  12.0   3.906  0.0198
##  SHAM transition - DREADD treatment     -3.269 0.807  11.9  -4.050  0.0157
##  SHAM transition - SHAM treatment        0.252 0.182 763.0   1.387  0.7349
##  DREADD treatment - SHAM treatment       3.521 0.805  11.8   4.375  0.0093
## 
## Degrees-of-freedom method: satterthwaite 
## P value adjustment: tukey method for comparing a family of 6 estimates

CamkII-hM3D(Gq) Hurst exponent

# H
p1 = make_timeplot(data4plot=data, x_var="time", y_var="H_basenorm", grp_var="id",
                   yLabel="H",title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_H_timeplot.pdf"), plot = p1)
p1

# run LME
res = run_lme(data = data, dv_var = "H_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## grp            13.02  13.020     1  11.05  7.1671   0.02145 *  
## condition     116.98  58.489     2 763.00 32.1966 3.767e-14 ***
## grp:condition 195.11  97.556     2 763.00 53.7024 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 3.766969e-14
res$aov_tab[3,6]
## [1] 1.511591e-22
res$pairwise_res
## $emmeans
##  grp    condition  emmean    SE   df lower.CL upper.CL
##  DREADD baseline    0.000 0.383 14.3   -0.820    0.820
##  SHAM   baseline    0.000 0.303 14.3   -0.649    0.649
##  DREADD transition -0.694 0.375 13.2   -1.503    0.116
##  SHAM   transition  0.441 0.297 13.2   -0.200    1.081
##  DREADD treatment  -2.145 0.370 12.5   -2.949   -1.342
##  SHAM   treatment   0.398 0.293 12.5   -0.237    1.033
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                             estimate    SE    df t.ratio p.value
##  DREADD baseline - SHAM baseline        0.0000 0.489  14.3   0.000  1.0000
##  DREADD baseline - DREADD transition    0.6936 0.206 763.0   3.369  0.0103
##  DREADD baseline - SHAM transition     -0.4405 0.485  13.9  -0.909  0.9378
##  DREADD baseline - DREADD treatment     2.1454 0.197 763.0  10.898  <.0001
##  DREADD baseline - SHAM treatment      -0.3980 0.482  13.6  -0.825  0.9577
##  SHAM baseline - DREADD transition      0.6936 0.482  13.6   1.438  0.7055
##  SHAM baseline - SHAM transition       -0.4405 0.163 763.0  -2.707  0.0751
##  SHAM baseline - DREADD treatment       2.1454 0.479  13.2   4.483  0.0062
##  SHAM baseline - SHAM treatment        -0.3980 0.156 763.0  -2.557  0.1094
##  DREADD transition - SHAM transition   -1.1341 0.478  13.2  -2.370  0.2344
##  DREADD transition - DREADD treatment   1.4519 0.181 763.0   8.029  <.0001
##  DREADD transition - SHAM treatment    -1.0915 0.476  12.9  -2.293  0.2638
##  SHAM transition - DREADD treatment     2.5860 0.475  12.8   5.448  0.0013
##  SHAM transition - SHAM treatment       0.0426 0.143 763.0   0.298  0.9997
##  DREADD treatment - SHAM treatment     -2.5434 0.472  12.5  -5.386  0.0015
## 
## Degrees-of-freedom method: satterthwaite 
## P value adjustment: tukey method for comparing a family of 6 estimates
p2 = make_scatterplot(data4plot = data2use,
                      x_var = "r_basenorm",
                      y_var = "H_basenorm",
                      grp_var = "grp",
                      cols2use = cols2use, 
                      xLabel = "Firing Rate",
                      yLabel = "H",
                      title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_H_scatterplot.pdf"), plot = p2)
p2

res = cor.test(data2use$r_basenorm, data2use$H_basenorm)
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm and data2use$H_basenorm
## t = -10.554, df = 323, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5830460 -0.4207607
## sample estimates:
##        cor 
## -0.5063734
res$p.value
## [1] 1.441197e-22
res = cor.test(data2use$r_basenorm[data2use$grp=="DREADD"], 
         data2use$H_basenorm[data2use$grp=="DREADD"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm[data2use$grp == "DREADD"] and data2use$H_basenorm[data2use$grp == "DREADD"]
## t = -3.1517, df = 123, p-value = 0.002039
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4284003 -0.1026790
## sample estimates:
##        cor 
## -0.2733575
res$p.value
## [1] 0.002038969
res = cor.test(data2use$r_basenorm[data2use$grp=="SHAM"], 
         data2use$H_basenorm[data2use$grp=="SHAM"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm[data2use$grp == "SHAM"] and data2use$H_basenorm[data2use$grp == "SHAM"]
## t = -3.3387, df = 198, p-value = 0.001005
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.35813380 -0.09517126
## sample estimates:
##        cor 
## -0.2308639
res$p.value
## [1] 0.001005481

CamkII-hM3D(Gq) 1/f slope

# 1/f slope
p3 = make_timeplot(data4plot=data, x_var="time", y_var="slope_basenorm", grp_var="id",
                   yLabel="1/f slope",title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_slope_timeplot.pdf"), plot = p3)
p3

# run LME
res = run_lme(data = data, dv_var = "slope_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
## grp           0.6579 0.65793     1  11.11  0.6340 0.4426
## condition     0.6974 0.34872     2 763.00  0.3360 0.7147
## grp:condition 4.0365 2.01827     2 763.00  1.9447 0.1437
res$aov_tab[2,6]
## [1] 0.7147202
res$aov_tab[3,6]
## [1] 0.1437322
res$pairwise_res
## $emmeans
##  grp    condition   emmean    SE   df lower.CL upper.CL
##  DREADD baseline    0.0000 0.201 19.9   -0.419    0.419
##  SHAM   baseline    0.0000 0.159 19.9   -0.331    0.331
##  DREADD transition  0.1994 0.192 16.7   -0.206    0.605
##  SHAM   transition -0.1811 0.152 16.7   -0.502    0.140
##  DREADD treatment   0.0171 0.186 14.8   -0.381    0.415
##  SHAM   treatment  -0.1305 0.147 14.8   -0.445    0.184
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                             estimate    SE    df t.ratio p.value
##  DREADD baseline - SHAM baseline        0.0000 0.256  19.9   0.000  1.0000
##  DREADD baseline - DREADD transition   -0.1994 0.156 763.0  -1.281  0.7954
##  DREADD baseline - SHAM transition      0.1811 0.252  18.6   0.720  0.9770
##  DREADD baseline - DREADD treatment    -0.0171 0.149 763.0  -0.115  1.0000
##  DREADD baseline - SHAM treatment       0.1305 0.249  17.9   0.524  0.9944
##  SHAM baseline - DREADD transition     -0.1994 0.249  17.9  -0.800  0.9637
##  SHAM baseline - SHAM transition        0.1811 0.123 763.0   1.472  0.6822
##  SHAM baseline - DREADD treatment      -0.0171 0.245  16.7  -0.070  1.0000
##  SHAM baseline - SHAM treatment         0.1305 0.118 763.0   1.110  0.8774
##  DREADD transition - SHAM transition    0.3805 0.245  16.7   1.555  0.6364
##  DREADD transition - DREADD treatment   0.1823 0.137 763.0   1.334  0.7662
##  DREADD transition - SHAM treatment     0.3299 0.242  15.9   1.363  0.7471
##  SHAM transition - DREADD treatment    -0.1982 0.240  15.5  -0.824  0.9584
##  SHAM transition - SHAM treatment      -0.0506 0.108 763.0  -0.468  0.9972
##  DREADD treatment - SHAM treatment      0.1476 0.238  14.8   0.621  0.9876
## 
## Degrees-of-freedom method: satterthwaite 
## P value adjustment: tukey method for comparing a family of 6 estimates
p4 = make_scatterplot(data4plot = data2use,
                      x_var = "r_basenorm",
                      y_var = "slope_basenorm",
                      grp_var = "grp",
                      cols2use = cols2use, 
                      xLabel = "Firing Rate",
                      yLabel = "1/f slope",
                      title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_slope_scatterplot.pdf"), plot = p4)
p4

res = cor.test(data2use$r_basenorm, data2use$slope_basenorm)
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm and data2use$slope_basenorm
## t = 1.3504, df = 323, p-value = 0.1778
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.03414132  0.18223593
## sample estimates:
##        cor 
## 0.07492923
res$p.value
## [1] 0.1778205
res = cor.test(data2use$r_basenorm[data2use$grp=="DREADD"], 
         data2use$slope_basenorm[data2use$grp=="DREADD"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm[data2use$grp == "DREADD"] and data2use$slope_basenorm[data2use$grp == "DREADD"]
## t = 2.3733, df = 123, p-value = 0.01918
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.03492914 0.37121966
## sample estimates:
##       cor 
## 0.2092531
res$p.value
## [1] 0.01918032
res = cor.test(data2use$r_basenorm[data2use$grp=="SHAM"], 
         data2use$slope_basenorm[data2use$grp=="SHAM"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm[data2use$grp == "SHAM"] and data2use$slope_basenorm[data2use$grp == "SHAM"]
## t = -0.078323, df = 198, p-value = 0.9377
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1441957  0.1332779
## sample estimates:
##          cor 
## -0.005566076
res$p.value
## [1] 0.9376503

CamkII-hM3D(Gq) Broadband power

# broadband power
p5 = make_timeplot(data4plot=data, x_var="time", y_var="isp_basenorm", grp_var="id",
                   yLabel="Power",title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_power_timeplot.pdf"), plot = p5)
p5

# run LME
res = run_lme(data = data, dv_var = "isp_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
##                Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## grp             4.529   4.529     1  11.1  3.5611   0.08557 .  
## condition     121.880  60.940     2 763.0 47.9183 < 2.2e-16 ***
## grp:condition  37.754  18.877     2 763.0 14.8435  4.74e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 2.490172e-20
res$aov_tab[3,6]
## [1] 4.740382e-07
res$pairwise_res
## $emmeans
##  grp    condition   emmean    SE   df lower.CL upper.CL
##  DREADD baseline    0.0000 0.236 18.5   -0.494   0.4940
##  SHAM   baseline    0.0000 0.186 18.5   -0.391   0.3906
##  DREADD transition -0.4196 0.226 15.8   -0.900   0.0609
##  SHAM   transition -0.0173 0.179 15.8   -0.397   0.3626
##  DREADD treatment  -1.4785 0.221 14.3   -1.951  -1.0059
##  SHAM   treatment  -0.3848 0.175 14.3   -0.758  -0.0112
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                             estimate    SE    df t.ratio p.value
##  DREADD baseline - SHAM baseline        0.0000 0.300  18.5   0.000  1.0000
##  DREADD baseline - DREADD transition    0.4196 0.172 763.0   2.436  0.1453
##  DREADD baseline - SHAM transition      0.0173 0.296  17.4   0.058  1.0000
##  DREADD baseline - DREADD treatment     1.4785 0.165 763.0   8.976  <.0001
##  DREADD baseline - SHAM treatment       0.3848 0.293  16.8   1.313  0.7745
##  SHAM baseline - DREADD transition      0.4196 0.293  16.8   1.431  0.7089
##  SHAM baseline - SHAM transition        0.0173 0.136 763.0   0.127  1.0000
##  SHAM baseline - DREADD treatment       1.4785 0.289  15.8   5.119  0.0012
##  SHAM baseline - SHAM treatment         0.3848 0.130 763.0   2.955  0.0378
##  DREADD transition - SHAM transition   -0.4023 0.289  15.8  -1.394  0.7301
##  DREADD transition - DREADD treatment   1.0589 0.151 763.0   6.999  <.0001
##  DREADD transition - SHAM treatment    -0.0348 0.286  15.2  -0.122  1.0000
##  SHAM transition - DREADD treatment     1.4612 0.284  14.8   5.142  0.0014
##  SHAM transition - SHAM treatment       0.3675 0.120 763.0   3.072  0.0266
##  DREADD treatment - SHAM treatment     -1.0937 0.281  14.3  -3.887  0.0161
## 
## Degrees-of-freedom method: satterthwaite 
## P value adjustment: tukey method for comparing a family of 6 estimates
p6 = make_scatterplot(data4plot = data2use,
                      x_var = "r_basenorm",
                      y_var = "isp_basenorm",
                      grp_var = "grp",
                      cols2use = cols2use, 
                      xLabel = "Firing Rate",
                      yLabel = "Power",
                      title2use="CamkII-hM3D(Gq)")
ggsave(file.path(plotdir,"camk_power_scatterplot.pdf"), plot = p6)
p6

res = cor.test(data2use$r_basenorm, data2use$isp_basenorm)
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm and data2use$isp_basenorm
## t = -3.8348, df = 323, p-value = 0.0001511
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3104191 -0.1022020
## sample estimates:
##        cor 
## -0.2086741
res$p.value
## [1] 0.0001511161
res = cor.test(data2use$r_basenorm[data2use$grp=="DREADD"], 
         data2use$isp_basenorm[data2use$grp=="DREADD"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm[data2use$grp == "DREADD"] and data2use$isp_basenorm[data2use$grp == "DREADD"]
## t = -0.6226, df = 123, p-value = 0.5347
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2293995  0.1207461
## sample estimates:
##         cor 
## -0.05604991
res$p.value
## [1] 0.5346984
res = cor.test(data2use$r_basenorm[data2use$grp=="SHAM"], 
         data2use$isp_basenorm[data2use$grp=="SHAM"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm[data2use$grp == "SHAM"] and data2use$isp_basenorm[data2use$grp == "SHAM"]
## t = 1.6355, df = 198, p-value = 0.1035
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.0236671  0.2501865
## sample estimates:
##      cor 
## 0.115453
res$p.value
## [1] 0.103532

hM4Di Firing Rate

data_file = here("experimental_data","hm4d_databasenorm.csv")
data = read.csv(data_file)
data$grp[data$grp=="exp"] = "DREADD"
data$grp[data$grp=="ctrl"] = "SHAM"
data$grp = factor(data$grp)
data2use = data %>% 
  filter(condition=="treatment")

# Firing Rate
p0 = make_timeplot(data4plot=data, x_var="time", y_var="r_basenorm", grp_var="id",
                   yLabel="Firing Rate",title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_FR_timeplot.pdf"), plot = p0)
p0

# run LME
res = run_lme(data = data, dv_var = "r_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
##                Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)    
## grp             4.959   4.959     1   8.01  2.5621 0.1481    
## condition     157.418  78.709     2 586.00 40.6677 <2e-16 ***
## grp:condition 244.156 122.078     2 586.00 63.0757 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 2.891814e-17
res$aov_tab[3,6]
## [1] 1.54989e-25
res$pairwise_res
## $emmeans
##  grp    condition  emmean    SE   df lower.CL upper.CL
##  DREADD baseline    0.000 0.692 8.69    -1.57    1.573
##  SHAM   baseline    0.000 0.692 8.69    -1.57    1.573
##  DREADD transition -1.777 0.687 8.46    -3.35   -0.208
##  SHAM   transition -0.337 0.687 8.46    -1.91    1.232
##  DREADD treatment  -2.836 0.684 8.32    -4.40   -1.268
##  SHAM   treatment   0.328 0.684 8.32    -1.24    1.895
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                             estimate    SE     df t.ratio p.value
##  DREADD baseline - SHAM baseline         0.000 0.978   8.69   0.000  1.0000
##  DREADD baseline - DREADD transition     1.777 0.213 586.00   8.363  <.0001
##  DREADD baseline - SHAM transition       0.337 0.975   8.57   0.346  0.9991
##  DREADD baseline - DREADD treatment      2.836 0.203 586.00  13.956  <.0001
##  DREADD baseline - SHAM treatment       -0.328 0.973   8.50  -0.337  0.9992
##  SHAM baseline - DREADD transition       1.777 0.975   8.57   1.823  0.4994
##  SHAM baseline - SHAM transition         0.337 0.213 586.00   1.587  0.6073
##  SHAM baseline - DREADD treatment        2.836 0.973   8.50   2.914  0.1286
##  SHAM baseline - SHAM treatment         -0.328 0.203 586.00  -1.612  0.5908
##  DREADD transition - SHAM transition    -1.440 0.972   8.46  -1.482  0.6833
##  DREADD transition - DREADD treatment    1.059 0.187 586.00   5.671  <.0001
##  DREADD transition - SHAM treatment     -2.105 0.970   8.39  -2.171  0.3394
##  SHAM transition - DREADD treatment      2.498 0.970   8.39   2.577  0.2037
##  SHAM transition - SHAM treatment       -0.665 0.187 586.00  -3.562  0.0053
##  DREADD treatment - SHAM treatment      -3.163 0.968   8.32  -3.269  0.0805
## 
## Degrees-of-freedom method: satterthwaite 
## P value adjustment: tukey method for comparing a family of 6 estimates

hM4Di Hurst exponent

# H
p1 = make_timeplot(data4plot=data, x_var="time", y_var="H_basenorm", grp_var="id",
                   yLabel="H",title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_H_timeplot.pdf"), plot = p1)
p1

# run LME
res = run_lme(data = data, dv_var = "H_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
##               Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)    
## grp             4.36   4.362     1   8.01   1.779  0.219    
## condition     320.59 160.294     2 586.00  65.376 <2e-16 ***
## grp:condition 196.79  98.394     2 586.00  40.130 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 2.348727e-26
res$aov_tab[3,6]
## [1] 4.637511e-17
res$pairwise_res
## $emmeans
##  grp    condition  emmean    SE   df lower.CL upper.CL
##  DREADD baseline    0.000 0.867 8.55   -1.977     1.98
##  SHAM   baseline    0.000 0.867 8.55   -1.977     1.98
##  DREADD transition  2.089 0.862 8.36    0.116     4.06
##  SHAM   transition  0.157 0.862 8.36   -1.816     2.13
##  DREADD treatment   3.296 0.859 8.25    1.325     5.27
##  SHAM   treatment   0.402 0.859 8.25   -1.569     2.37
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                             estimate    SE     df t.ratio p.value
##  DREADD baseline - SHAM baseline         0.000 1.226   8.55   0.000  1.0000
##  DREADD baseline - DREADD transition    -2.089 0.239 586.00  -8.733  <.0001
##  DREADD baseline - SHAM transition      -0.157 1.223   8.46  -0.128  1.0000
##  DREADD baseline - DREADD treatment     -3.296 0.229 586.00 -14.410  <.0001
##  DREADD baseline - SHAM treatment       -0.402 1.221   8.40  -0.329  0.9993
##  SHAM baseline - DREADD transition      -2.089 1.223   8.46  -1.708  0.5601
##  SHAM baseline - SHAM transition        -0.157 0.239 586.00  -0.655  0.9866
##  SHAM baseline - DREADD treatment       -3.296 1.221   8.40  -2.700  0.1729
##  SHAM baseline - SHAM treatment         -0.402 0.229 586.00  -1.757  0.4945
##  DREADD transition - SHAM transition     1.932 1.219   8.36   1.585  0.6275
##  DREADD transition - DREADD treatment   -1.207 0.210 586.00  -5.745  <.0001
##  DREADD transition - SHAM treatment      1.687 1.217   8.31   1.386  0.7346
##  SHAM transition - DREADD treatment     -3.139 1.217   8.31  -2.579  0.2039
##  SHAM transition - SHAM treatment       -0.245 0.210 586.00  -1.167  0.8524
##  DREADD treatment - SHAM treatment       2.894 1.215   8.25   2.381  0.2632
## 
## Degrees-of-freedom method: satterthwaite 
## P value adjustment: tukey method for comparing a family of 6 estimates
p2 = make_scatterplot(data4plot = data2use,
                      x_var = "r_basenorm",
                      y_var = "H_basenorm",
                      grp_var = "grp",
                      cols2use = cols2use, 
                      xLabel = "Firing Rate",
                      yLabel = "H",
                      title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_H_scatterplot.pdf"), plot = p2)
p2

res = cor.test(data2use$r_basenorm, data2use$H_basenorm)
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm and data2use$H_basenorm
## t = -13.91, df = 248, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7264211 -0.5860897
## sample estimates:
##        cor 
## -0.6620185
res$p.value
## [1] 6.649365e-33
res = cor.test(data2use$r_basenorm[data2use$grp=="DREADD"], 
         data2use$H_basenorm[data2use$grp=="DREADD"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm[data2use$grp == "DREADD"] and data2use$H_basenorm[data2use$grp == "DREADD"]
## t = -14.142, df = 123, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.8456471 -0.7093022
## sample estimates:
##       cor 
## -0.786895
res$p.value
## [1] 1.481409e-27
res = cor.test(data2use$r_basenorm[data2use$grp=="SHAM"], 
         data2use$H_basenorm[data2use$grp=="SHAM"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm[data2use$grp == "SHAM"] and data2use$H_basenorm[data2use$grp == "SHAM"]
## t = -9.99, df = 123, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.7560316 -0.5594230
## sample estimates:
##        cor 
## -0.6692811
res$p.value
## [1] 1.440809e-17

hM4Di 1/f slope

# 1/f slope
p3 = make_timeplot(data4plot=data, x_var="time", y_var="slope_basenorm", grp_var="id",
                   yLabel="1/f slope",title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_slope_timeplot.pdf"), plot = p3)
p3

# run LME
res = run_lme(data = data, dv_var = "slope_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
##                Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## grp            3.9455  3.9455     1   8.06  2.9629    0.1232    
## condition      0.3827  0.1913     2 586.00  0.1437    0.8662    
## grp:condition 28.3783 14.1891     2 586.00 10.6556 2.848e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 0.8661809
res$aov_tab[3,6]
## [1] 2.847621e-05
res$pairwise_res
## $emmeans
##  grp    condition  emmean    SE    df lower.CL upper.CL
##  DREADD baseline    0.000 0.257 12.54  -0.5575    0.557
##  SHAM   baseline    0.000 0.257 12.54  -0.5575    0.557
##  DREADD transition -0.358 0.248 10.91  -0.9047    0.189
##  SHAM   transition  0.229 0.248 10.91  -0.3181    0.776
##  DREADD treatment  -0.570 0.243  9.99  -1.1114   -0.029
##  SHAM   treatment   0.524 0.243  9.99  -0.0171    1.065
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                             estimate    SE     df t.ratio p.value
##  DREADD baseline - SHAM baseline         0.000 0.364  12.54   0.000  1.0000
##  DREADD baseline - DREADD transition     0.358 0.176 586.00   2.029  0.3271
##  DREADD baseline - SHAM transition      -0.229 0.357  11.71  -0.640  0.9853
##  DREADD baseline - DREADD treatment      0.570 0.169 586.00   3.383  0.0099
##  DREADD baseline - SHAM treatment       -0.524 0.354  11.23  -1.482  0.6814
##  SHAM baseline - DREADD transition       0.358 0.357  11.71   1.001  0.9088
##  SHAM baseline - SHAM transition        -0.229 0.176 586.00  -1.299  0.7859
##  SHAM baseline - DREADD treatment        0.570 0.354  11.23   1.612  0.6074
##  SHAM baseline - SHAM treatment         -0.524 0.169 586.00  -3.109  0.0240
##  DREADD transition - SHAM transition    -0.587 0.351  10.91  -1.670  0.5749
##  DREADD transition - DREADD treatment    0.213 0.155 586.00   1.373  0.7434
##  DREADD transition - SHAM treatment     -0.882 0.347  10.45  -2.539  0.1969
##  SHAM transition - DREADD treatment      0.799 0.347  10.45   2.301  0.2740
##  SHAM transition - SHAM treatment       -0.295 0.155 586.00  -1.906  0.3991
##  DREADD treatment - SHAM treatment      -1.094 0.343   9.99  -3.186  0.0774
## 
## Degrees-of-freedom method: satterthwaite 
## P value adjustment: tukey method for comparing a family of 6 estimates
p4 = make_scatterplot(data4plot = data2use,
                      x_var = "r_basenorm",
                      y_var = "slope_basenorm",
                      grp_var = "grp",
                      cols2use = cols2use, 
                      xLabel = "Firing Rate",
                      yLabel = "1/f slope",
                      title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_slope_scatterplot.pdf"), plot = p4)
p4

res = cor.test(data2use$r_basenorm, data2use$slope_basenorm)
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm and data2use$slope_basenorm
## t = 6.6842, df = 248, p-value = 1.523e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2802265 0.4909770
## sample estimates:
##       cor 
## 0.3907097
res$p.value
## [1] 1.523323e-10
res = cor.test(data2use$r_basenorm[data2use$grp=="DREADD"], 
         data2use$slope_basenorm[data2use$grp=="DREADD"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm[data2use$grp == "DREADD"] and data2use$slope_basenorm[data2use$grp == "DREADD"]
## t = 4.898, df = 123, p-value = 2.98e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2458295 0.5412082
## sample estimates:
##       cor 
## 0.4039967
res$p.value
## [1] 2.980145e-06
res = cor.test(data2use$r_basenorm[data2use$grp=="SHAM"], 
         data2use$slope_basenorm[data2use$grp=="SHAM"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm[data2use$grp == "SHAM"] and data2use$slope_basenorm[data2use$grp == "SHAM"]
## t = 2.1629, df = 123, p-value = 0.03249
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.01635512 0.35508475
## sample estimates:
##       cor 
## 0.1914129
res$p.value
## [1] 0.03248631

hM4Di Broadband power

# broadband power
p5 = make_timeplot(data4plot=data, x_var="time", y_var="isp_basenorm", grp_var="id",
                   yLabel="Power",title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_power_timeplot.pdf"), plot = p5)
p5

# run LME
res = run_lme(data = data, dv_var = "isp_basenorm")
res$aov_tab
## Type III Analysis of Variance Table with Satterthwaite's method
##                Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## grp             25.67   25.67     1   8.02  1.4412    0.2642    
## condition     2073.29 1036.64     2 586.00 58.1915 < 2.2e-16 ***
## grp:condition  409.89  204.95     2 586.00 11.5045 1.257e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
res$aov_tab[2,6]
## [1] 8.867731e-24
res$aov_tab[3,6]
## [1] 1.256753e-05
res$pairwise_res
## $emmeans
##  grp    condition  emmean   SE   df lower.CL upper.CL
##  DREADD baseline    0.000 1.53 9.36   -3.448     3.45
##  SHAM   baseline    0.000 1.53 9.36   -3.448     3.45
##  DREADD transition  4.000 1.51 8.90    0.570     7.43
##  SHAM   transition  0.498 1.51 8.90   -2.933     3.93
##  DREADD treatment   6.643 1.50 8.62    3.223    10.06
##  SHAM   treatment   2.634 1.50 8.62   -0.787     6.05
## 
## Degrees-of-freedom method: satterthwaite 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                             estimate    SE     df t.ratio p.value
##  DREADD baseline - SHAM baseline         0.000 2.168   9.36   0.000  1.0000
##  DREADD baseline - DREADD transition    -4.000 0.645 586.00  -6.205  <.0001
##  DREADD baseline - SHAM transition      -0.498 2.155   9.13  -0.231  0.9999
##  DREADD baseline - DREADD treatment     -6.643 0.616 586.00 -10.776  <.0001
##  DREADD baseline - SHAM treatment       -2.634 2.146   8.99  -1.227  0.8143
##  SHAM baseline - DREADD transition      -4.000 2.155   9.13  -1.857  0.4795
##  SHAM baseline - SHAM transition        -0.498 0.645 586.00  -0.772  0.9722
##  SHAM baseline - DREADD treatment       -6.643 2.146   8.99  -3.095  0.0961
##  SHAM baseline - SHAM treatment         -2.634 0.616 586.00  -4.272  0.0003
##  DREADD transition - SHAM transition     3.503 2.141   8.90   1.636  0.5981
##  DREADD transition - DREADD treatment   -2.643 0.566 586.00  -4.667  0.0001
##  DREADD transition - SHAM treatment      1.367 2.132   8.76   0.641  0.9843
##  SHAM transition - DREADD treatment     -6.145 2.132   8.76  -2.882  0.1320
##  SHAM transition - SHAM treatment       -2.136 0.566 586.00  -3.772  0.0024
##  DREADD treatment - SHAM treatment       4.009 2.124   8.62   1.888  0.4663
## 
## Degrees-of-freedom method: satterthwaite 
## P value adjustment: tukey method for comparing a family of 6 estimates
p6 = make_scatterplot(data4plot = data2use,
                      x_var = "r_basenorm",
                      y_var = "isp_basenorm",
                      grp_var = "grp",
                      cols2use = cols2use, 
                      xLabel = "Firing Rate",
                      yLabel = "Power",
                      title2use="hM4Di")
ggsave(file.path(plotdir,"hM4Di_power_scatterplot.pdf"), plot = p6)
p6

res = cor.test(data2use$r_basenorm, data2use$isp_basenorm)
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm and data2use$isp_basenorm
## t = -8.017, df = 248, p-value = 4.275e-14
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5469561 -0.3492670
## sample estimates:
##        cor 
## -0.4536751
res$p.value
## [1] 4.275038e-14
res = cor.test(data2use$r_basenorm[data2use$grp=="DREADD"], 
         data2use$isp_basenorm[data2use$grp=="DREADD"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm[data2use$grp == "DREADD"] and data2use$isp_basenorm[data2use$grp == "DREADD"]
## t = -6.0163, df = 123, p-value = 1.892e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.6020279 -0.3287518
## sample estimates:
##        cor 
## -0.4768312
res$p.value
## [1] 1.892143e-08
res = cor.test(data2use$r_basenorm[data2use$grp=="SHAM"], 
         data2use$isp_basenorm[data2use$grp=="SHAM"])
res
## 
##  Pearson's product-moment correlation
## 
## data:  data2use$r_basenorm[data2use$grp == "SHAM"] and data2use$isp_basenorm[data2use$grp == "SHAM"]
## t = -4.4494, df = 123, p-value = 1.904e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.5143217 -0.2105001
## sample estimates:
##        cor 
## -0.3723438
res$p.value
## [1] 1.903682e-05

Make figures of power spectrum

freq_data = read.csv(here("experimental_data","freq.csv"))
freq_names = freq_data$colnames2use

# CamkII=-hM3D(Gq)
data_file = here("experimental_data","LFPspectra_xsub_camk.csv")
data = read.csv(data_file)
mean_data = data %>% filter(condition=="treatment")

# n1 = length(unique(mean_data$id[mean_data$grp=="SHAM"]))
n1 = sum(mean_data$grp=="SHAM")
data1 = data.frame(mean_spectrum_bn = colMeans(mean_data[mean_data$grp=="SHAM",freq_names]), grp = "SHAM", freq = freq_data$freq, se_spectrum_bn = (apply(mean_data[mean_data$grp=="SHAM",freq_names], 2, sd) / sqrt(n1)))

# n2 = length(unique(mean_data$id[mean_data$grp=="DREADD"]))
n2 = sum(mean_data$grp=="DREADD")
data2 = data.frame(mean_spectrum_bn = colMeans(mean_data[mean_data$grp=="DREADD",freq_names]), grp = "DREADD", freq = freq_data$freq, se_spectrum_bn = (apply(mean_data[mean_data$grp=="DREADD",freq_names], 2, sd) / sqrt(n2)))

data2plot = rbind(data1,data2)
data2plot$logFreq = log10(data2plot$freq)

camk_pow_spec = ggplot(data = data2plot, 
                       aes(x = logFreq,
                           y = mean_spectrum_bn,
                           colour=grp,
                           fill=grp)) + 
  geom_ribbon(aes(ymin = mean_spectrum_bn - se_spectrum_bn,
                  ymax = mean_spectrum_bn + se_spectrum_bn), 
              alpha = 0.2) + 
  geom_line() + 
  geom_hline(yintercept=0) +
  xlab("Frequency (Hz)") + 
  ylab("Baseline Normalized Power") + 
  scale_x_continuous(labels = c(0.1, 1, 10, 100)) +
  ggtitle("CamkII-hM3D(Gq)") + easy_center_title() +
  easy_add_legend_title("Group") + easy_adjust_legend(to = "center") +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, size=fontSize),
        strip.text.y = element_text(size=fontSize),
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize-2),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text = element_text(size = fontSize))




# hM4Di
data_file = here("experimental_data","LFPspectra_xsub_hm4d.csv")
data = read.csv(data_file)
mean_data = data %>% filter(condition=="treatment")

# n1 = length(unique(mean_data$id[mean_data$grp=="SHAM"]))
n1 = sum(mean_data$grp=="SHAM")
data1 = data.frame(mean_spectrum_bn = colMeans(mean_data[mean_data$grp=="SHAM",freq_names]), grp = "SHAM", freq = freq_data$freq, se_spectrum_bn = (apply(mean_data[mean_data$grp=="SHAM",freq_names], 2, sd) / sqrt(n1)))

# n2 = length(unique(mean_data$id[mean_data$grp=="DREADD"]))
n2 = sum(mean_data$grp=="DREADD")
data2 = data.frame(mean_spectrum_bn = colMeans(mean_data[mean_data$grp=="DREADD",freq_names]), grp = "DREADD", freq = freq_data$freq, se_spectrum_bn = (apply(mean_data[mean_data$grp=="DREADD",freq_names], 2, sd) / sqrt(n2)))

data2plot = rbind(data1,data2)
data2plot$logFreq = log10(data2plot$freq)

hm4di_pow_spec = ggplot(data = data2plot, 
                       aes(x = logFreq,
                           y = mean_spectrum_bn,
                           colour=grp,
                           fill=grp)) + 
  geom_ribbon(aes(ymin = mean_spectrum_bn - se_spectrum_bn,
                  ymax = mean_spectrum_bn + se_spectrum_bn), 
              alpha = 0.2) + 
  geom_line() + 
  geom_hline(yintercept=0) +
  xlab("Frequency (Hz)") + 
  ylab("Baseline Normalized Power") + 
  scale_x_continuous(labels = c(0.1, 1, 10, 100)) +
  ggtitle("hM4Di") + easy_center_title() +
  easy_add_legend_title("Group") + easy_adjust_legend(to = "center") +
  theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5, size=fontSize),
        strip.text.y = element_text(size=fontSize),
        axis.text.x = element_text(size=fontSize),
        axis.text.y = element_text(size=fontSize-2),
        axis.title.x = element_text(size=fontSize),
        axis.title.y = element_text(size=fontSize),
        strip.text = element_text(size = fontSize))

p_final = camk_pow_spec / hm4di_pow_spec + plot_layout(guides = "collect")
ggsave(filename = file.path(plotdir, "dreadd_pow_spec.pdf"), 
       plot = p_final,
       width = 8,height=8)
p_final